Bose condensates in a harmonic trap near the critical temperature 
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The mean-field properties of finite-temperature Bose-Einstein gases confined in spherically sym- 
metric harmonic traps are surveyed numerically. The solutions of the Gross-Pitaevskii (GP) and 
Hartree-Fock-Bogoliubov (HFB) equations for the condensate and low-lying quasiparticle excitations 
are calculated self-consistently using the discrete variable representation, while the most high-lying 
states are obtained with a local density approximation. Consistency of the theory for temperatures 
through the Bose condensation point T c requires that the thermodynamic chemical potential dif- 
fer from the eigenvalue of the GP equation; the appropriate modifications lead to results that are 
continuous as a function of the particle interactions. The HFB equations are made gapless either 
by invoking the Popov approximation or by renormalizing the particle interactions. The latter ap- 
proach effectively reduces the strength of the effective scattering length a sc , increases the number of 
condensate atoms at each temperature and raises the value of T c relative to the Popov approxima- 
tion. The renormalization effect increases approximately with the log of the atom number, and is 
most pronounced at temperatures near T c . Comparisons with the results of quantum Monte Carlo 
calculations and various local density approximations are presented, and experimental consequences 
are discussed. 
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I. INTRODUCTION 

Since the first observations of Bose-Einstein condensa- 
tion (BEC) in dilute alkali metal atom gases ex- 
perimental developments have posed many new tests for 
many-body theory, even though weakly interacting Bose 
gases have long been used as a textbook paradigm ||,^| . 
Numerous theoretical approaches have been employed in 
order to obtain accurate results for both the ground-state 
and non-equilibrium properties of the trapped Bose sys- 
tems However, there have been notable differences 
between theoretical results and experimental data on the 
excitation frequencies near the transition temperature 
T c ]9|-|I1|. This problem has inspired the introduction 
of a renormalized effective atom-atom interaction [ pd| . 
Recently developed theoretical approaches [ [l3|]i!| ] that 
incorporate the dynamics of the noncondensate density 
but without a renormalized interaction have resulted in 
excitation frequencies in closer agreement with experi- 
ment. Nevertheless, the unresolved issues for Bose sys- 
tems near T c has provided a motivation for us to examine 
further the theoretical and numerical methods for model- 
ing confined Bose gases near T c . We have numerically im- 
plemented the most plausible and tractable equilibrium 
mean-field theories in order to systematically survey var- 
ious properties of these systems. 

In this work, we follow the standard mean-field the- 
ory |l5| , with certain modifications described in detail 
below. The nonlinear Gross-Pitaevskii (GP) equation, 
which includes interactions between the condensate and 
the thermal atoms, is solved for a static condensate con- 
taining N c atoms. The eigenvalue of the GP equation, 



ft, is usually identified with the thermodynamic chemical 
potential fi. The linear response of the system is repre- 
sented by the Hartree-Fock-Bogoliubov (HFB) equations, 
which yield the quasiparticle energies and amplitudes. 
These in turn determine the number of noncondensed 
atoms Nt as well as various coherence terms (thermo- 
dynamic averages over two or more Bose field opera- 
tors). The GP and HFB equations are iterated to self- 
consistency at a given temperature T, subject to a fixed 
total number of atoms in the system N = N c + Nt- 
As emphasized by Griffin [jl5| , the coherence terms yield 
an excitation spectrum that is not gapless: the lowest- 
energy mode of the HFB equations has finite energy and 
does not coincide with the solution of the GP equation. 
The HFB-Popov (HFBP) approximation, which neglects 
these terms, has been quite successful in describing the 
properties of the trapped Bose gases, but is not well- 
grounded theoretically, and fails to yield accurate pre- 
dictions for the low-lying excitations at high tempera- 
tures |9|,[T^| . In this work, we explore a recently proposed 
extension of the HFBP theory that incorporates the co- 
herence terms in a gapless manner Jl6| , pT| ; in addition, we 
modify the commonly used identification of the chemical 
potential with the eigenvalue of the GP equation. 

The identification of the chemical potential with the 
eigenvalue of the GP equation is incorrect in general. 
In the grand canonical ensemble, the chemical poten- 
tial is defined as — dE/dN, corresponding to the en- 
ergy cost E of adding a particle to the entire system, not 
only to the condensate. For a dilute, weakly interacting 
Bose gas at T = 0, for which the population of noncon- 
densed states (the depletion) is negligible, the identifica- 
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tion p, = fj, is justified. At finite temperatures, however, 
the assumption yields results that are discontinuous as a 
function of the s-wave scattering length a sc . To a bet- 
ter approximation, we find that the chemical potential 
at finite temperatures is given by the eigenvalue of the 
GP equation plus a term that varies inversely with the 
number of condensate atoms. The resulting equations 
provide an improved description of these finite systems, 
yielding observables that are both continuous with a sc 
and similar to those obtained using path integral Monte 
Carlo techniques JL7[ . 

It is presently unclear to what extent many-body ef- 
fects beyond the mean-field approximation modify the ef- 
fective interactions am ong Bose-condensed atoms in har- 
monic traps jy],|l6|[ij|jl£] . For the homogeneous Bose 
gas, it is now established both from rcnormalization 
group |^0| and perturbation pl[ theories that the many- 
body T-matrix, or effective s-wave scattering length a, 
goes to zero at T c . The low-energy, long wavelength limit 
of the many-body T-matrix has been shown to be closely 
related to the coherence term tut |l6| , |l9| ; this 'anomalous 
average' represents two-particle correlations and is the 
Bose analogue of the superconducting order parameter 
in interacting Fermi systems. Renormalizing the interac- 
tion using the mj yields a gapless HFB theory without 
having to invoke the Popov approximation |Tl[| , but it re- 
mains uncertain whether the prescription is appropriate 
for nonuniform systems. The implications of this theory 
for trapped Bose condensates are explored numerically 
below, and the results are compared to those obtained 
with the Popov approximation and path-integral Monte 
Carlo methods. 

In view of these somewhat conflicting results and un- 
resolved issues, there is strong motivatation for contin- 
ued development of numerical methods in order to imple- 
ment various models and obtain quantitative predictions 
for comparison with experiment. Quantum Monte Carlo 
methods jlT],^],^ are able to provide accurate results for 
certain observable quantities. The computational proce- 
dure is lengthy, however, and is not demonstrably able to 
yield excitation frequencies since it typically applies only 
to equilibrium configurations. Local density approxima- 
tions (LDA) are much simpler to apply, but the standard 
forms fail near T c and are questionable when the density 
is so small that the local collision rate is insufficient to 
establish local thermodynamic equilibrium. On the other 
hand, widely used basis set techniques are generally un- 
able to represent the large numbers of atoms in excited 
states at high-temperatures. Recently, Reidl et al. J|4|] 
have used (for 2,000 Rb atoms at T = 0.5T C ) a hybrid 
method in which a sum over discrete quasiparticle states 
at low energies is supplemented by an integral over an 
energy-dependent LDA above some cutoff energy. The 
interactions of these two subensembles with each other 
are expressed by mean-field potentials that represent the 
effect of background atoms. In the present work, the low- 
lying states are obtained by solving the HFB equations 
using the discrete variable representation (DVR) 1 25 2?]] 



and the cutoff energy is raised until the results converge 
to within a stated tolerance. The techniques employed 
have enabled the investigation of trapped Bose gases at 
finite-temperatures containing a larger number of atoms 
than in previous calculations that we are aware of. As a 
result, the approach of these systems to the local ther- 
modynamic equilibrium and to the hydrodynamic limit 
can be explored. 

In Section II A, we outline the GP and HFB equations. 



We disc uss th e ch emic al potential and gaple ss the ories in 
Sections |IIB| and |ll Cj , respectively. Section II D reviews 
LDA methods both as alternative approaches for compar- 
ison purposes, and the com ple mentary use for the most 
energetic atoms. In Section III, we discuss our numerical 
methods and iteration procedures. Section |v] presents 
results for Bose atoms in a spherically symmetric har- 
monic trap as a function of the scaled s-wave scattering 
length, total number of atoms, and temperature. 



II. THEORETICAL FRAMEWORK 

A. Thermal sums over quasiparticle states 

The derivation of mean-field equations for a weakly 
interacting, dilute Bose gas has been described in de- 
tail elsewhere |2^-|3^Jl^] . The question of the chemical 
potential for T > for thermal sums of quasiparticle 
states deserves more thorough discussion, however, and 
we modify the standard procedure. In addition, follow- 
ing discussions by the Burnett et al. p6|,|i"l 19 1, we treat 



the anomalous (coherence) terms m-r in a manner that 
produces a 'gapless' theory. 

Following the standard approach, we decompose the 
Bose field operator into a c-number for the condensate, 
plus an operator representing its fluctuations. The full 
many-body Hamiltonian is approximated using mean- 
field theory, becoming explicitly number-nonconserving. 
The grand canonical ensemble is used, and thus the chem- 
ical potential, //, and temperature, T, are the sole fixed 
quantities. The generalized Gross-Pitaevskii (GP) equa- 
tion for the condensate and coupled Bogoliubov equa- 
tions for the excited quasiparticle states are then solved. 
For a finite number atoms in a harmonic potential, how- 
ever, the standard approach yields values for the mean 
condensate number N c that are discontinuous as a func- 
tion of interaction strength a sc . In our approach, the 
eigenvalue of the GP equation, jl, is determined by the 
mean number of atoms in the condensate N c . In con- 
trast, the chemical potential, fi, is adjusted so that the 
mean total number of atoms is the desired value. A sim- 
ple relationship is found connecting jl, [i, and N c , which 
is adapted from the ideal Bose gas case. 

The Hamiltonian for an interacting Bose gas in a trap 
in the grand canonical ensemble is 



H-/J.N 



dr 



2M 



V 2 + V a 
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(1) 



where the field operator ip(r) satisfies [ip(ri), ^(r?)] = 
8(yi — r 2 ). The pseudopotential atom-atom interaction 
has been chosen to be V{r\ — r 2 ) = g5(ri — where 
the coupling constant g = 4Trh 2 a sc /M is written in terms 
of the scattering length a sc and mass M. The harmonic 
potential is V ex t — hMoi^r 2 with trapping frequency ujo 
is assumed to be isotropic. 

The Hamiltonian may be rewritten as 



H - fiN = H - [IN + (/x - fj)N 
= K+ (/i - n) N, 



(2) 



where, as mentioned above, the Lagrange multiplier fl 
is related to the number of atoms in the condensate. 
In the following, we choose to diagonalize the opera- 
tor K = H — p,N rather than the original H — /j,N; 
both choices must lead to the same excitation spectrum, 
though with a temperature-dependent shift of the vac- 
uum for quasiparticle excitations. In order to make fur- 
ther progress, the Bose field operator if) — $ + <f> is now 
decomposed into a c-number for the condensate, and 
(f>(r), which annihilates a thermal atom at r. The con- 
densate density is defined by n c — |$| 2 , and the number 
of condensate atoms is N c = f dr\&(r)\ 2 . The noncon- 
densate (thermal) density tit and anomalous (coherence) 
terms my and mj are ]l5| ] 

n T = m T = rh T = (<^ t t ), (3) 



where the brackets indicate a thermal average in the 
grand canonical ensemble, discussed in more detail be- 
low. The mean field approximation is used to reduce the 
third and fourth order terms to, respectively, first and 
second order in <j>, cj>' ! so that the Hamiltonian K can be 
diagonalized, following the procedure normally used for 

Excluding the possibility of aggregate motion and vor- 
tices [E9[, <3? may be taken to be real. The first order 
terms (plus third order terms in mean- field approxima- 
tion) in K vanish if the equation for the condensate is 
taken to be the generalized GP equation: 



2M 



V 2 + V ext + g[n c + 2n T + fh T ] 



$ = (4) 



Note that /2 is the eigenvalue of the GP equation. The 
part of K that is zeroth order in the excited orbitals is a 
c-number 



dr$(r) 



2M 



fJ 



$(r)| 2 $(r). 



(5) 



The terms in K that are second order in cj> are (in the 
mean-field approximation) diagonalized by the canonical 
transformation 



4>(r) = ^[ Ui (r)ay + v*(r)a]] 

3 



(0) 



such that [oijjCit] = Sij. The operator K is diagonal to 

second order in <j> if the quasiparticle amplitudes Uj(r) 
and Vj(r) are solutions of the Bogoliubov coupled equa- 
tions 



Cuj(r) + Q(r)vj(r) = e u {v) 
Cvj{r) + Q(r) Uj (r) = -e^-(r), 



(7) 



where t = K + V cxt - £ + 2 .9™( r ); Q = .9K(r) + m T {r)], 
the total density is n(r) = n c (r) + nr(r). For 'gapless' 
theories, discussed further below, the j = 'Goldstone 
mode,' has the property eo = 0, so that Uq(t) = — «o( r ) = 
<f>(r). Thus on the Ej energy scale, the condensate has 
zero energy, and defines the vacuum for quasiparticle ex- 
citations. 

After the substitution, ip = $ + <p, the number opera- 
tor N — J dvi/ji (r)^>(r) contains terms, such as J dr$<f>, 
that do not conserve particle number. The Bogoliubov 
transformation (||) and coupled equations (Q) introduce 
a quasiparticle basis such that terms ctjCtj and ctj&j 
are eliminated, so quasiparticle number is conserved p8[ . 
The diagonalized Hamiltonian explicitly does not con- 
serve particle number, however; the operator K in the 
quasiparticle basis does not commute with the excited 
particle number operator (jMj), which has contributions 
from cvcv and aa terms. In the grand canonical en- 
semble, only T and /i are precisely defined, and all ob- 
servables must be defined in terms of thermal averages. 
Each occupation number, including the condensate num- 
ber, fluctuates about its mean value 



{Nj) = (a' a), j = 0,1, 



(8) 



where the explicit definition of the average (O) is yet 
undefined. Similarly, both the eigenvalue fl of the GP 
equation ([|) and the total energy (E) fluctuate about 
their mean values. 

Inserting the transformation (^) into Eqs. (|^) and in- 
troducing the identification given by Eq. (|8|) , the normal 
and anomalous densities become 

Mr) = £{<#i>[|«i(r)| a + Mr)| 2 ] + Mr)| 2 } i (9) 

3=1 



m T (r) = ^^(r)vt(r)[2(^) + 1]. (10) 

3=1 

The standard normalization J dr[\uj(r)\ 2 — 1 1;^ (r ) | ] = 1, 
yields 



dr[\u j (r)f + \v j (r)\ 2 ] = l + 2V j , 



(11) 
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where Vj — J dr\vj(r)\ 2 . The quantities Vj are related 
to the T = depletion, which is Yj- The relation 

between the total atom number and the quasiparticle oc- 
cupation numbers is therefore 



(12) 



(N) =N C + N T = (N ) + J dvn T (v) 
= N c + Y l [(N j )(l + 2V j ) + V j ], 



where the average number of atoms is written in terms 
of a contribution from the condensate and noncondensate 
(excitations). 

The thermal average of the diagonalized Hamiltonian 
then becomes 

(H - fxN) =K +Y, ejdNj) - Vj) + (N) 
= E c - nN c + { + (£-/*)(! + Wj)] 

3=1 

+V i {fi-n-e j )}, (13) 

where E c = Kq + jlN c is the total ground state, or con- 
densate, energy. 



Eq. (|l4|), deviations of Nj> from (Nj') for j' ^ j will not 
greatly modify the spectrum of the quasiparticle states. 
If this is so, a reasonable approximation is to replace 
(Nj) by Nj when estimating the mean occupation num- 
bers from the grand partition function. If the dependence 
of ej and Nj on Nj/ (j ^ j') is also neglected, then fi can 
be factored, and we obtain 



(Nj) 



Eat, Nj exp {-/? [ej + (£ - //)(! + 2Vj)] Nj} 
J2 Nj exp {-f3 [ej + (/i - M )(l + 2Vj)] Nj} 



exp{/3 [e J + (ft-^(l + 2V J )}}-l' 



Vj. (16) 



In order to obtain the result on the second line of 
Eq. (|l~6|), the population-dependences of the GP eigen- 
value jl and the quasiparticle energies tj are ignored. At 
sufficiently low temperatures, the tj for trapped Bose- 
condensates are relatively insensitive to the value of N c 
and the temperature; indeed, in the Thomas-Fermi (TF) 
limit, valid for large condensates, the excitation frequen- 
cies at zero temperature are independent of N c . 

Neglecting the factors Vj , and shifting the energy scale 
so that Ej = ej + fM, one recovers the more conventional 
expression 



B. Occupation factors and the chemical potential 

In the Bogoliubov approach, the ensemble is considered 
to be the sum of a condensate plus non-interacting quasi- 
particles. The mean occupation numbers of the quasipar- 
ticle states are to be determined from the grand partition 
function, 



n = Tr{exp [-/3(H - fiN)]}, 
through the standard identities |^,|| 

/An 1 /<91nQ\ /ainfi 



(14) 



(15) 



Unfortunately, while the diagonalized Hamiltonian is 
written in terms of non- interacting single-quasiparticle 
energies, the expressions (|l^) and (jlj) involve the ther- 
mal averages of particle occupation that we are now seek- 
ing to determine. Furthermore, the factorization that one 
makes for an ideal Bose gas is invalid for a gas of interact- 
ing Bose atoms because the quasiparticle energies depend 
on the occupation numbers, as well as the reverse. Thus, 
rigorously, these occupation factors should be calculated 
self-consistently, along with Eqs. (§) and (fjj), since they 
depend on as well as determine the quasiparticle eigenval- 
ues |n| . To do so analytically would be a truly daunting 
task. We make several simplying assumptions in order 
to obtain results, but we emphasize that these questions 
merit further study. 

In reality, the probabilities (N) will be peaked at the 
most probable values, as discussed below for the conden- 
sate. Therefore, when evaluating the sum over Nj in 



exp[p(Ej - a*)] - 1 ' 



(17) 



From this expression 

@ for 3 = 0, with E = ji, one 
finds that the chemical potential /i and the eigenvalue of 
the GP equation jl are related by the expression 



In 1 



1 



N c >0. 



(18) 



For T = this gives the usual definition jl = [i, but for 
T > there is a correction to \i that increases as N c 
decreases. While this additional term will not be correct 
at high temperatures where the condensate is strongly 
depleted, it will be shown below that results obtained 
with this procedure are continuous functions of a sc at all 
temperatures, while with /i = jl they are not. 

It is difficult to go beyond the above approximations, 
but we will suggest possible avenues to proceed in future 
work. The major effect omitted is the dependence of the 
quasiparticle energies, Ej (including Eq = jl) on N c . One 
can first consider the condensate term itself. We assume, 
for the moment, that factorization of SI ( |l4| ) is valid, and 
write 



(19) 



In the Thomas-Fermi approximation (kinetic energy in 
the GP equation neglected), one obtains for a spherical 
condensate [p|: 



Mtf 



1 n5N c a s N 2/5 

2 V a 



JN, 



2/5 



(20) 
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where the harmonic oscillator length is oq — y/h/Muj. 
The following relations follow in the same approximation: 



dE c 
dN r _ 



(21) 



Then from Eq. (p_9[), neglecting Qt, one can obtain 
the mean value of the condensate occupation from 
(N c ) = En c n cP(Nc)/En c P(N c ), where P(N C ) = 

exp[/3(/ziV c — j^Nc )]■ We have verified numerically for 
typical values of (3 and /i that the mean value is extremely 
close to the most probable value N c for which P(N C ) 
is maximum. Furthermore, an expansion of the expo- 
nent in the above expression for P(N C ) yields a value 
for the variance of N e , interpreted as the value of a in 
such that P(N C ± a) — (l/e)P(N c ). In the grand canon- 
ical ensemble at zero temperature, therefore, one obtains 
(SN C ) = y/(5/P"/)N? /W , so that the fractional width of 

the occupation number distribution decreases as N c 7 / 10 . 
This may be compared with the result of Giorgini et al., 
derived from excited state occupation numbers for the 
canonical ensemble, (SN C ) ~ (T/T c )iV c 2/3 ||. Either re- 
sult confirms that the fluctuations in 7V C are relatively 
small for large N c . One should next consider how the 
dependence of Qt would effect (N c ) and {SN C ). This is 
left for future work. 

The dependence of the quasiparticle states on the oc- 
cupation factors reflects the extensive nature of this fi- 
nite interacting system; that is, adding a particle to the 
many-body system alters both the number and character 
of the accessible states. This behavior is similar | fj3|]34| 
to that of a finite gas of non-interacting particles obeying 
fractional exclusion statistics []35|| , which obey a statistics 
intermediate between that of bosons and fermions. The 
parameter representing the statistics has been identified 
with the strength of the delta-function potential for an 
interacting trapped Bose gas in two dimensions |54|. In- 
deed, our expression ([l8]) for the thermodynamic chemi- 
cal potential is similar to that found for a non-interacting 
fractional-statistics gas at finite temperature [^3|,^4| ■ We 
hope to pursue these issues more fully in future work. 



C. Gapless approximations 

We return to the conditions for 'gaplessness.' The 
GP (|J) and Bogoliubov (Q) equations together comprise 
the 'Hartree-Fock-Bogoliubov' (HFB) approximation for 
a dilute interacting Bose gas. In this case, one does not 
obtain e = 0, and the theory is said to be not gapless 
(the term has been taken from the homogeneous situ- 
ation). In the Popov approximation, gaplessness is en- 
sured by neglecting the coherence terms tut and tjit, but 
the justification for such an approximation is question- 
able [£|. 

In order to convert HFB into a gapless theory and still 
retain the anomalous averages, Burnett et al. [ |16|Jllf| have 



recently proposed an alternative treatment in which the 
coupling functions for the condensate g c (r) and excited 
states g e (r) absorb the pairing correlations, and thereby 
take on a spatial dependence. Eq. (Q) becomes 



{K + V cxt + g c n c + 2sf e n T }$ = /i4>, 



(22) 



and similarly, £ and Q appearing in Eqs. (|7|) become 

C = K + V CKt - M + ^g c n c + 2g e nT; Q = g c nc- (23) 

In the proposed gapless theories, labeled Gl, and G2, the 
coupling constants are chosen to be: 



where 



{9c: 9e\ 



3i W 



{31:3}, Gl 
{31; Si}, G2 



mr(r) 
n c (r) 



(24) 



(25) 



The renormalized coupling g\ replaces the two-body T- 
matrix associated with binary atomic collisions, which 
is the scattering length a sc in vacuo, by the zero mo- 
mentum and energy limit of the homogeneous many- 
body T-matrix (T^|ll],[ji| . In the Gl approximation, only 
the condensate-condensate and condensate-excited are 
dressed, while G2 is motivated by the expectation that 
all particle interactions should be similar. Renormaliza- 
tion of the coupling has the additional advantage of re- 
moving the ultraviolet divergence in my resulting from 
high-energy quasiparticle contributions of Eq. (|o|) in the 
T-matrix approximation. In nonuniform systems such 
as the trapped Bose gases, however, the value of g\ (r) 
can diverge in regions near the condensate surface where 
the condensate density vanishes more rapidly than the 
anomalous average. In practice, this divergence may be 
eliminated by setting g\(r) = g [1 + mT{v)/ (n c (r) + 5)], 
where S rj 1CP 2 . While the results, described in detail 
below, are found not to depend strongly on the choice of 
5, its existence underlines a deficiency in the theory in 
its present form. The consequences of the Gl approxi- 
mation are not explored in this work. In the following, 
the notation g(r) will be used in place of <?i(r) and in 
distinction to g, which is unrenormalized. 



D. Local Density Approximation 

In local density approximation (LDA) schemes, the 
condensate density is assumed to be varying sufficiently 
slowly that the population of excited states is determined 
entirely by the local potential and temperature. The 
thermal density may then be treated locally as if the 
interacting Bose gas were homogeneous. We will discuss 
three basic LDA schemes, and several variants. 
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In the semiclassical approximation to the GP and HFB 
equations (24|]3(|, the thermal atom quasiparticle ampli- 
tudes in the Bogoliubov equations (0) become local func- 
tions u(p, r) and u(p, r). With the Popov approximation, 
one obtains the coupled algebraic equations 



£(p,r) gn c (v) \ ( it(p,r) 
-gn c (r) -C(p, r) J\ v(p, r) 



e(p,r) 



u(p,r) 
u(p,r) 



(26) 



where £(p, r) = p 2 /2m + V e ^t( r ) — A + 2.gn(r). With 
the condition u(p,r) 2 — u(p, r) 2 = 1, the local excita- 
tion energies may be immediately obtained, e(p, r) = 
(£(p, r) 2 — g 2 n 2 (r) 2 ) 1 / 2 , and have the well known lin- 
ear dispersion. The noncondcnsate density from Eq. (j^) 
may then be easily found p4[ : 



n T (r) 



d 3 p 



£(P,r) 
e (P, r ) 



;n(p,r)) 



e(£( P ,r) 2 -.g 2 n 2 (r)) 



where 



(n(p,r)) 



1 



exp[/3(e(p,r) + /i — jLt)] — 1 ' 



(27) 



(28) 



such that the theta function is unity when the argument 
is positive, and zero otherwise. These equations define 
the Hartree-Fock Bogoliubov Popov LDA, which we will 
refer to as the "BPLDA." For G2 calculations, one ob- 
tains the "BGLDA" by the substituion g — > ^(r) every- 
where. Then one needs 



mr(r) 



d 2 p 

— -r^(p,r)u(p,r)[2(n(p,r)) + f] 



;(r)n c (r) 



d 3 p 1 



(2tt) 3 2e 
9(/:(p,r) 2 - 5 2 nc(r)) 



[2(n(p,r)) + l] 



(29) 



The integral is not formally convergent, however. Since 
the anomalous averages appear only in the context of the 
Gl and G2 approximations, where the formal ultraviolet 
divergence is eliminated, we may safely neglect the +1 
term following the 2(n(p, r)). 

The semiclassical HFBP approximation exhibits a gap- 
less excitation spectrum only if the condensate is also 
treated within the LDA, which implies the TF density: 



n c {r) = 



jj. - Vext(r) ~ 2gn T (r) 

g 

e[/i-V 0Xt (r)-2 5 n T (r)], 



(30) 



The TF approximation is valid in the limit of large 
N c , where the energy contribution from the mean-field 
(Hartree) potential exceeds that of the kinetic energy. 
For this reason, Eq. ( |30| ) is not expected to be a good 
representation of the condensate density close to the tran- 
sition temperature. 



In the regime of small condensate numbers, therefore, 
it becomes more important to solve the equations for 
the condensate and excitations exactly in order to obtain 
the low-lying discrete states, as described in the previous 
section. In this work, we use the exact GP and HFB 
equations, but the sum over discrete states is combined 
with an energy integral over high-lying states usin g L DA 
functions in the manner described by Reidl et al. P4fl : 



n T (r) 



E 



nj(r)Q(e c - ej) + den T (e,r), 

J e„ 



(31) 



where rij(r) is the jth term of Eq. (||), e c is a low-energy 
cutoff, and, in the above notation, ny(e, r) has the form 



n T (e,r) 



,3/2 



2 3 / 2 7T 

x [C - T4xt 



2(n(p,r)) + l-- 



Ji - 2gn] 



1/2 



(32) 



A similar equation applies to the anomalous average mx- 
This latter hybrid procedure is referred to below as the 
Discrete Quasiparticle Sum (DQS) approximation, an ab- 
breviation for 'discrete Hartree-Fock-Bogoliubov quasi- 
particle sum.' Either a Popov or G2 approximation may 
be made within the DQS, and these are referred to below 
as DQSP and DQSG, respectively. 

A simpler LDA may be formulated by treating the lo- 
cal excitations within the Hartree-Fock approximation, 
which ignores the linear dispersion at low energies. The 
condensate density may again be obtained within the TF 
approximation using Eq. (|30|). The thermal density is 

given by n T (r) = J ■^ s {n(p,r)), where (n(p,r)> is de- 
fined in Eq. ( p8|) but with e(p,r) = C(p, r). Integration 
over the momenta readily yields 



n T (r) 



1 



ff3/2 



-/3(Voxt(r)+2 ff n(r)- /1 ) 



(33) 



where the thermal de Broglie wavelength is At = 
(2tt^ 2 /mkT) 1 ' 2 and g a {z) = Y,f =1 z s /j a . As usual, the 
chemical potential /i is determined by the condition that 
the total atomic number, N = J dr[n c (r) +iit(v)]. With 
the TF expression (|30| ) for the condensate, the argument 
of the c/3/2 function in Eq. ( |33"| ) is always less than unity. 
If an 'exact' solution for the condensate is used (i.e. ob- 
tained by solving the GP equation), the results are gener- 
ally improved, but as noted below and in Ref . , there 
is then a range of temperatures T % T c for which the 
g%/2 function given in Eq. ( |33"| ) diverges, since its argu- 
ment can become greater than unity. 

An even simpler form of the LDA has been formu- 
lated |3^^^1 in which the effect of interactions on the ex- 
cited states is completely ignored. Assuming a TF form 
for the ground state, this LDA consists of the parametri- 
cally coupled equations (in view of the other approxima- 
tions here, in these equations we ignore the distinction 
between /2 and [i): 
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n c (r) 



M - V"cxt(r) 
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nr(r) = ^3-53/2 

Arp 



9[/t-F ext (r)] 

-/3|Voxt(r)-Ar 



(34) 
(35) 



In this approximation, the interaction enters only via the 
chemical potential in the TF equation, which is a function 
of a sc and condensate number. For a spherical conden- 
sate, /EtxF = ^(15N c a BC /ao) 2 ^ 5 huio, where ao = ^/K/Mujq 
is the bare oscillator length. 

It is shown in Ref. |37[] that a low order expansion of 
Eq. (|35l) yields the following expression for N /N: 



Nc = (t_\ _ cm 

N {t°J \(3) 



T 
T9 



N 



2/5 



(36) 



where C( n ) is the Riemann zeta function, r\ = 
fi TF /k B T° re iC(3) 1/3 (15A^ 1/6 a sc /a ) 2/5 , and the criti- 
cal temperature for N ideal Bose atoms in a harmonic 
trap is given by |3§|]3l| 

k B T°/huj = 0.9405A 1/3 - 0.6842 + 0.50A~ 1/3 (37) 

Equation (|36j) is solved iteratively for N c /N. 

E. Ideal Bose gas 

Some of the plots given below contain results for ideal 
non-interacting Bose atoms (a sc = 0) in a harmonic trap. 
The results given for N c were obtained from sums over 
the occupation numbers as given in Eq. (|J), with dj = 
2£j + l, Ej = hu{t+2rij+Z/2). The chemical potential \i 
was adjusted to satisfy the condition N = Ylj=a(Nj}- An 
alternative expression can be obtained from the density 
distribution given by Chou et al. Eol : 
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l-Zl 
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e 2 ^)(tanh(/3V2)]~ 3/2 -l}, (38) 



where Zi 



_ „/3(/x-3/2) 



This expression requires fewer 



terms than the aforementioned procedure, and gives iden- 
tical results for temperatures up to about 0.9T C . 



III. COMPUTATIONAL TECHNIQUES 

With a spherically symmetric trapping potential, all 
observables may be decomposed into functions of radius 
r and spherical harmonics yj n {9, 4>). The GP and Bogoli- 
ubov equations then become one-dimensional in r; the 
ground state is assumed to have (£, m) = (0,0), while 
the excitations obtained using the Bogoliubov equations 
are 2£+l degenerate. Both equations are solved using the 



discrete variable representation (DVR), a computation- 
ally efficient approach for the trapped interacting Bose 
gases that has been recently described in detail |2^] . 

We have used two variants of the DVR approach: an 
equidistant mesh array derived from sine functions as dis- 
cussed by Colbert and Miller (^6), and a mesh based on 
Gaussian quadrature, using the zeros of associated La- 
guerre polynomials L% L if) , where Nl is the order of the 
the quadrature and a — 2 for a spherical condensate |27| . 
The latter DVR has the advantage of having a fine mesh 
for small r where the condensate density is non-zero, and 
a more coarse mesh at larger distances where the ther- 
mal distribution varies smoothly. Though the condensate 
and excited orbitals are computed on the physical grid, 
the matrix elements of the operators are represented by 
Laguerre polynomials up to the order defining the Gauss 
quadrature Nl, which in the present calculations range 
from 1,000 to 2,800; matrix elements of the kinetic en- 
ergy are computed from expressions given in Ref. [p7| . 
Increasing the value of Nl increases the accuracy of high- 
lying states, allowing for a larger cutoff energy e c at which 
the discrete sums are terminated, and a smaller number 
of atoms in the LDA integrals. Since high-order polyno- 
mials extend far beyond values of R max ^ 50ao relevant 
to trapped condensates, the number of spatial grid points 
required can be limited to just N g ~ 200 for all values of 
Ni,. 

Implementation of the above mean-field theory re- 
quires a stable and efficient iteration procedure to solve 
the GP and Bogoliubov equations for a given total num- 
ber of atoms iV and temperature T . In our approach, 
the functions n c (r) — <f> 2 (r) , tit (t) , and nriT{r) are cal- 
culated self-consistently using Eqs. (||) and (|7|)-(fi0|), sup- 
plemented by the LDA expressions for states above the 
cutoff e c , for fixed N c and T; the chemical potential \x is 
determined by Eq. (|l8|). Because this iteration procedure 
is especially delicate near T c , yet is crucial for the results 
presented, we give a few more details. 

We emphasize that the convergence criterion must con- 
sider the spatial distribution functions n c (r) and nr(r) 
rather than simply the aggregate values, N c and Nt- The 
iterative procedure can be decomposed into three sepa- 
rate levels of self-consistency, subject to the minimization 
of the 'Error': 



Error 



dr[|n° ut (r) - <(r)| + \n^(r) - n£(r)|]. 



(39) 



The 'in' and 'out' functions are the input and output 
of the combined GP and HFB equations plus the high 
energy LDA integral. Normally, the Error diminishes 
(though not necessarily monotonically) through level 1 
iterations, in which the output functions are fed back into 
the GP, HFB and high energy LDA equations. In this 
level, the condensate number N c is held constant while 
the condensate density (normalized to unity) is allowed 
to vary. When the Error reaches some predetermined 
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tolerance, level 2 iterations begin and N c is adjusted to 
approach the condition that N c + Nt = N. The first 
level 2 adjustment from the converged level one itera- 
tions is based on a simple proportionality between N 
and N c . Subsequent level 2 adjustments are based on 
a linear relation between A?" c and N, where the param- 
eters are obtained from the last two level 2 iterations. 
After N c + Nt has converged to N to the desired toler- 
ance, level 3 iterations proceed, in which iteration levels 
1 and 2 are repeated with successively larger number of 
Laguerre functions Nl and mesh points N g . These three 
levels of iteration typically achieve accuracies for the con- 
densate number N c of a few atoms. While this accuracy 
is beyond what is accessible to current experiments, it 
permits the comparison of different theoretical models. 
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1. Convergence of the self-consistency procedure, for 
10 5 , a BC /a = 0.0072, t BC = 53, and a Laguerre DVR 



tion of atoms in the LDA integral, Fi a t- (c) Cut-off energy, e c 
specifying the upper limit of the discrete quasiparticle sum. 
(d) Order of the Laguerre polynomial, Nl. (e) Condensate 
number N c . (f) Error, defined by Eq. ([ii^), showing conver- 
gence up to each change of N c or N g , and ultimately conver- 
gence to the condition that N c + Nt = N. 

The iteration procedure is illustrated in Fig. [I], which 
tracks a calculation for 2 • I0 5 atoms and scaled temper- 
ature t = k B T/Tiuj = 53 (from Eq. @, t° « 54.3), using 



the Laguerre DVR basis. After more than 50 iterations, 
N c converged from the initial estimate of 109 to the fi- 
nal value of 149 atoms (Fig. le). Each adjustment of 
N c (level 2) or Nt (level 3) results in a jump in the error 
(Fig. If), which then converges again. In this calculation, 
Nl increased from 1300 to 2100 (Fig. Id), corresponding 
to an increase of mesh points (up to i? max = 42) from 
149 to 190 (Fig. la), an increase in e c from 102^0 to 
lA4hujo (Fig. lc), and a decrease in the fraction of total 
number of atoms in the LDA integral from 57% to 40% 
(Fig. lb). 

The fraction of atoms in the LDA integral is negligi- 
ble only for calculations at low temperatures with small 
N. Since T c rises as ~ 0.947V 1 / 3 , the required number of 
thermal states rises with N for calculations near T c , and 
inevitably the LDA integration must include a larger frac- 
tion of atoms. For N = 2 • 10 4 , 2 ■ 10 5 m and 10 6 , at most 
9%, 38%, and 74% of the atoms were in the integral at 
temperatures in the vicinity of T c . Correspondingly, the 
mesh size N g required to ensure convergence increased 
from 140 to 210 for N between 10 3 and 10 6 . The reason 
N g does not increase more rapidly with N is that the 
LDA approximation improves with the total number of 
atoms. 

It should be pointed out that for large values of N, 
the iteration procedure could exhibit instabilities when 
the temperature approached T c . For N > 10 5 , we found 
that there often appeared to be (at least) two semi-stable 
regions when N c < 5, 000, between which the calculation 
tended to fluctuate. In order to ensure the solution re- 
mained in the more stable state, small temperature in- 
crements At = 0.2 were used. 



IV. THERMAL AVERAGES 

A. Condensate fraction 

In several of the plots to follow, results are presented 
for a series of values of a sc /aQ. For comparison with 
current experiments, we note that the scattering lengths 
a sc for 87 Rb, 23 Na, and 7 Li are approximately given by 
llOas, 52a,B, and — 27.3as, respectively, where as ~ 
5.292 ■ 10~ n m is the Bohr radius. Thus, if one takes u = 
(uixUjyUJz) 1 / 3 , then for the recent MIT experiments [12] 
on 23 Na, v = lu/2it = 96.4 Hz, the JILA experiments | ] 



mesh, (a) Number of points in the DVR mesh, N g . (b) Frac- g; 



ive v = 182.5 Hz, and the Rice experiments |3j give v = 
144.6 Hz, corresponding to a sc /a = 0.00129, 0.00729, 
and —0.00046, respectively. 

Figure [2] illustrates the consequences of setting the 
eigenvalue of the GP equation jl equal to the chemical 
potential //, as discussed in Section II A. With this as- 



sumption (here used in conjunction with the Popov ap- 
proximation, rriT = 0), N c goes to zero abruptly with 
T when the population in excited states reaches the to- 
tal number of atoms N = 5, 000. By contrast, results for 
a, r = 0, obtained as described in Sec. HE, have a smooth 
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tail at high temperature. Thus, in the limit a sc — > 0, the 
results near T c exhibit a discontinuity with respect to the 
ideal gas results. 
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FIG. 2. If fi — ft, from the HFBP discrete quasiparticle 
sum (DQS) near T c there is a discontinuity in the N c vs. T 
function with respect to a sc - The figures show TV C vs. T for 
TV = 5,000 atoms, for several values of a sc /ao- Even in the 
limit of small a BC /ao, N c goes to zero abruptly with T for the 
self-consistent solution, while for the ideal Bose gas (a sc = 0) , 
N C (T) has a smooth tail. 





24.8 



25.2 



24.4 
t sc = kT/tico 

FIG. 3. The chemical potential in units of hui relative to 
the harmonic oscillator zero-point energy, fi/hu — 3/2, vs. 
T for various values of a sc /ao. (a) Shows the full range of 
temperatures up to T c , while (b) shows a limited range near 
T c . 

Figures H and || show results obtained from calculations 
in which the chemical potential is as given in Eq. (18). 
The smooth variation of the chemical potential, Fig7]3|, 
through T c is reflected in all relevant properties of the 
system, including the number of condensate atoms and 
excitation frequencies. When a sc > 0, the chemical po- 
tential evolves continuously from positive to negative val- 
ues, relative to the harmonic oscillator zero point energy 
fftw, as the temperature increases. Since \i increases with 
the interaction strength, the value at which the chemi- 



cal potential passes through zero increases with a sc even 
though T c decreases. In addition, Fig. || shows that for 
a sc < 0, /J, < 3hui/2 everywhere, with maximum values 
at temperatures T ~ T c . 

Figure |] shows the number of condensate atoms as a 
function of temperature for TV = 1, 000 and 20, 000 for a 
range of interaction strengths a sc /ao, calculated within 
the DQS formalism. The condensate population near T c 
is evidently a continuous function of both the scattering 
length and temperature. 



200 




FIG. 4. When fi differs from fi according to Eq. (|18[), the 
Nc vs. T function from the discrete quasiparticle sum behaves 
smoothly with respect to a sc . Shown are the results for (a) 
TV = 1,000 and (b) TV = 20,000 atoms. The critical temper- 
ature for a ac = 0, defined as the maximum of d 2 N c /dT 2 , is 
indicated with an arrow. For a sc < the maximum value of 
TV C is limited due to the instability of the condensate. In (b) , 
open circles denote results obtained with the G2 approach. 

The plots shown in Fig. |J, especially for 20, 000 atoms, 
show that the G2 renormalization procedure results in 
a significantly higher value of TV C , relative to that ob- 
tained within the Popov approximation, for the larger 
values of a sc . Furthermore, the difference between the 
G2 and Popov results becomes more pronounced as a sc 
increases. This behavior is consistent with expectation 
because G2 produces a weakening of the atom-atom in- 
teraction. The use of the occupation factors Jig) rather 







than (g) also increases the value of N c by a few atoms at 
high temperatures, but the effect is much smaller than 
what results from the use of G2 theory. 

For a sc < 0, the N c values reach a maximum when 
the calculation becomes numerically unstable jll] B , 
reflecting the physical instability of the cloud towards 
spatial collapse. The maximum N c values depend on 
a sc , as shown by the termination of the curves for these 
cases. For T = 0, the maximum value is given by 
A™ ax = 0.573a /a sc @. This critical number is known 
to decrease when T > due to the presence of thermal 
atoms [^2[^3). In these plots, the maximum N c is 80% to 
57% of the value calculated for T = 0, confirming that 
the thermal cloud significantly decreases the stability of 
the condensate for a sc < 0. 



temperature values. 

It is remarkable that the values for N c from BLDA cal- 
culations agreed with the corresponding DQS results to 
better than 0.4% of N in every case for which results were 
obtained. Even for HFLDA and SILDA, the differences 
with DQS results are less than the fractional error in 
current experiments. Thus these comparisons show that 
relatively simple LDA expressions are useful for obtain- 
ing the condensate fraction as a function of temperature. 
It is only in the region near T c and above, where the con- 
densate number becomes small, that our LDA methods 
failed. 



B. Comparison with LDA and QMC 

It is interesting to explore how our finite temperature 
results compare with those obtained by other methods. 
Local density approximations are much simpler to imple- 
ment numerically than the full self-consistent HFB equa- 
tions and their variants. The opposite is true of Monte 
Carlo calculations, but these do not invoke the mean-field 
approximation so yield results for equilibrium configura- 
tions that are essentially exact. 

Fig. | for N = 2 • 10 4 compares N c values from the 
Popov and G2 quasiparticle sums (DQSP and DQSG) 
with severak LDA methods. Our Hartree-Fock LDA 
(HFLDA) solves the GP equation for the condensate 
n c (r), iterated to self-consistency using Eq. ( |33| ) for the 
thermal distribution nr(r). We found it most efficient 
to start at low temperature, in order to obtain good ini- 
tial estimates of nr(r) at successively higher values of 
T. No solution could be found for N c /N < 0.035 due 
to the failure of the HFLDA, as discussed above and in 
Ref. @. 

The 'semi-ideal' LDA (SILDA) J37| omits the n T (r) 
term in the TF expressions for the condensate (30) and 
for the total density ny(r) in Eq. ([331) . This results in the 
simple expressions ( pf ) which are related solely through 
the chemical potential. Iterative solution of these equa- 
tions yields results that are close to the other functions 
plotted in Fig. 0. The actual nr(r) distribution calcu- 
lated with this approach exhibits a sharp peak at edge of 
the condensate due to the discontinuity at the Thomas- 
Fermi condensate radius. 

The inset of Fig. || shows that the Hartree-Fock Bogoli- 
ubov LDA methods, BPLDA and BGLDA, agree most 
closely with the hybrid method, DQSP and DQSG, re- 
spectively. The two BLDA methods employ a TF con- 
densate, and thus the nr(r) functions exhibit a small 
spike at the edge of the condensate, which has a cusp. As 
with the HFLDA, the calculations required iteration to 
self-consistency, which was facilitated when initial values 
were obtained by extrapolation from results from lower 
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FIG. 5. Comparison of values for N c /N from quasiparti- 
cle sums with the Popov and G2 approximations, as com- 
pared with HFLDA and SILDA for N = 20, 000 atoms and 
tSsc/fflo = 0.0072. BLDA results are too close to distinguish on 
this scale. On an expanded scale, the inset gives differences 
between indicated LDA and DQS methods. 

The Quantum Monte Carlo (QMC) approach uses the 
exact Hamiltonian with a hard-sphere atom-atom inter- 
action. Based on extensive numerical experience with 
4 He p2] |, QMC should be most useful for the calculation 
of equilibrium quantities, such as the condensate fraction. 
Holzmann el al. jl7j have provided benchmark QMC cal- 
culations for the case of 10 4 Bose atoms confined in a 
spherical trap, with a sc /ao = 0.0043. Table | shows com- 
parisons between our results and those of QMC p7 45 
for the condensate number as a function of temperature. 
The DQSP, DQSG, and QMC values differ by up to 1.2% 
of the total atom number N. It is notable that at higher 
temperatures, N c falls off less quickly using HFBP and 
G2 than QMC. This may be due in part to the fact that 
the relationship between N c and /x in Eq. (|l^) is not 
entirely correct at higher temperatures, as discussed in 
Section |II B| , and may resemble ideal gas statistics too 
closely. Presumably the many-body effects that neces- 
sitate the renormalization of the atom-atom interaction 
are already included in the QMC procedure, in which 
case results with G2 should be closer than Popov to the 
QMC. Indeed, for t = k B T/Tiu < 17, the Popov results 
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lie below QMC, while the G2 numbers are higher and 
closer to QMC. Above a scaled temperature t = 18, how- 
ever, the G2 results rise above QMC values. 

C. Critical Temperature vs. a sc 

Figures ^ and || show that large values of a sc have the 
effect of flattening the curve of condensate number as a 
function of temperature, as is already apparent in the 
plots of Giorgini et al. pq ]. If these curves are fit to a 
function N c /N = 1 — (T/T c ) a , one obtains values for a 
as low as 1.4, compared with the ideal gas value of 3. An- 
other parameter to characterize the effect of atom-atom 
interactions is the shift of the critical temperature from 
the ideal Bose gas case. For the homogeneous Bose gas, 
where it is uniquely defined as the point at which N c 
goes to zero, this shift has been the subject of intense 
discussion recently |46|| . For atoms in a harmonic poten- 
tial, as is especially clear in Fig. [|, this point is not sharp 
(indeed, the number of condensate atoms is finite at all 
temperatures in a mesoscopic system). Definitions of T c 
that have been proposed include the point at which the 
density at the origin reaches the critical density for a ho- 
mogeneous gas pTJ], the maximum of the specific heat, 
and the maximum of the temperature derivative of the 
specific heat |39) . Since such energy- weighted properties 
pose additional problems for numerical calculations of 
thermal averages, T c is determined here as the maximum 
of the function d 2 N c /dT 2 . The inflection point of the N c 
vs. T function, or zero of d 2 N c /dT 2 deviated from Eq. 
( |37| ) by a significantly larger amount. 

Figureo shows T c values extracted from the data used 
in Fig. |P For comparison, the ideal gas data are an- 
alyzed in a similar manner, yielding values of T c that 
are close to, but not identical with, those obtained using 
Eq. (|37j). Figures ||a and ||b correspond to 1,000 and 
20, 000 atoms, respectively. The inset in Fig. ||b shows 
how the transition temperature is determined from the 
data in a typical case, by making use of the three func- 
tions N C (T), dN c (T)/dT, and d 2 N c /dT 2 . Since both the 
condensate number and its temperature derivative are 
nearly straight lines, accurate calculation of the second 
derivative requires accurate numeric values of these func- 
tions. 

A semiclassical analysis by Giorgini et al. [H indicates 
that the transition temperature should decrease linearly 
from the ideal gas value with increasing particle interac- 
tions. The results of the DQS-Popov calculations confirm 
this general scaling; furthermore, as the number of atoms 
increases, the observed shift in the critical temperature 
5T C matches the semiclassical expression more closely at 
larger a sc /ao- In contrast, with the DQS-G2 approach 
5T C shows significant deviations from linear scaling for 
small N, and these become more pronounced as the num- 
ber of atoms increases. For N = 2 ■ 10 4 , the shifts are sig- 
nificantly less than the semiclassical values for the larger 



values of a sc /ao considered. 
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FIG. 6. Values for the critical temperature, T c , defined 
as the maximum of d 2 N c /dT 2 . Results are shown for (a) 
N = 1,000 and (b) N = 20,000 atoms, for the DQSP and 
DQSG approaches. The solid line represents the semiclassical 
prediction T c = T° — AT C , where T° is the transition tem- 
perature in the non- interacting limit. The inset in (b) shows 
the N C (T), N c (T)/dT, and d 2 N c /dT 2 functions from which 
T c is determined for the cases a sc /ao = 0.0048, the last with 
a spline fit. 

D. Renormalization of the atom-atom interaction 

As indicated in Fig. [| above, the G2 renormaliza- 
tion yields values for N c /N that reflect the weakening 
of the atom-atom repulsion; at any given temperature, 
the number of atoms in the condensate increases rela- 
tive to the value obtained using the Popov approxima- 
tion. Perhaps more interesting is the spatial variation 
of the effective interaction in the harmonic trap. The 
renormalization is governed by the local value of tut rel- 
ative to n c . In general, |mr| increases with the number 
of noncondensed atoms ut since more terms enter the 
sum (|l0|); however, m-r vanishes when n c = 0, since the 
'quasihole' amplitude Vi = 0. In general, therefore, one 
might expect the local renormalized interaction to reach 
a minimum at some temperature. For a uniform Bose 
gas, this minimum occurs at exactly the transition tem- 
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perature, and corresponds to a vanishing of the effective 
scattering length 

In Fig. ^ we compare the condensate and thermal den- 
sities with the spatial variations of the anomalous av- 
erage and the effective particle interactions for the case 
of N = 20, 000 and a sc /ao = 0.0072 for various tem- 
peratures. (As noted above, this would correspond to a 
relatively tight trap for 23 Na.) For these plots, 6 = 0.01 
in Eq. (^||) . There is a slight dependence of the results 
on 6, since much smaller values S ~ 10~ 4 lead to a small 
bump in the mxir) function at the very edge of the con- 
densate. The dependence on this arbitrary parameter 
indicates an ambiguity in the theory; however, the inte- 
grated numbers are not significantly altered by the choice 
of 5, since the errors are incurred in regions of very small 
condensate density. 
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FIG. 7. The functions n c (r), nr(r), mr(r), and g(r)/g is 
shown for N = 20,000 atoms, a sc /ao = 0.0072 over a wide 
range of temperatures. |mT(r)| is largest at the edge of the 
condensate and increases with T up to T c . 

The manner in which g(r)/g attains a minimum in r 
is shown in Fig. || for the particular case of N = 2 • 10 5 . 
The global minimum occurs at a temperature close to 
T c , defined above. Following this procedure, we con- 
sider the <7min(j)/g functions for various values of N for 
a sc /ao — 0.0072, which are displayed in Fig. ||. Though 
we have increased N without changing the trap fre- 
quency, the approach to the thermodynamic limit is be- 
ginning to emerge. The minimum for each N is found 
to always occur very close to the calculated transition 
temperature, and its value decreases approximately with 
log(A) over the range of N considered. For N = 10 6 , we 
obtain g m in(r)/g ~ 0.2. It should be noted that although 
the fraction of total atoms in the LDA integral increased 
to approximately | for N = 10 6 near T c , the high-energy 
LDA contribution to was in every case less than 2%, 
and typically an order of magnitude less than this value. 




FIG. 8. Variation of the renormalization factor, g(r)/g 
with temperature near T c for TV = 200,000, and 
dsc/ao = 0.0072. The range of the minimum decreases as the 
condensate shrinks with T, while the minimum value contin- 
ues to decrease up to a point, and then increases. 

It should be emphasized that the G2 renormalization 
employed in the present calculations is derived for a uni- 
form Bose gas, and should best represent large conden- 
sate densities or low temperatures where the LDA is most 
applicable. While the LDA is bound to fail for T — > T c , 
the regime where it loses validity will become smaller 
with increasing A, and should approach the critical re- 
gion where perturbation theory itself breaks down. It 
would be preferable to define the renormalization of the 
particle interactions in terms of the full many-body T- 
matrix in a trap, and we hope to pursue this issue in fu- 
ture work. The G2 approach as formulated above, how- 
ever, should properly describe the effects of two-body 
correlations for large trapped condensates at low to in- 
termediate temperatures. Thus, the strong reduction in 
the effective interaction strength over much of the con- 
densate, indicated by the G2 theory, could have signifi- 
cant experimental consequences. The predictions for the 
excitation frequencies are discussed further below. 

E. Excitation frequencies 

The quasiparticle eigenvalues correspond to excitation 
frequencies, but it remains unclear what relationship ex- 
ists between these values and experimentally observed 
resonances of the trapped gas at finite temperatures when 
the potential of a harmonic trap is perturbed periodi- 
cally. In all mean-field calculations such as those pre- 
sented here, the linear response equations assume that 
the thermal density is fixed, while in experiments it would 
also be perturbed. For this reason, the dipole excitation 
frequency obtained within mean-field theories will gener- 
ally not satisfy the generalized Kohn theorem [jl9| , which 
states that there is a mode in which the entire ensemble 
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oscillates at the bare trap frequency. Calculations explic- 
itly including the dynamics of both n c and m [ |13|]l4|l are 
found to be consistent with the Kohn theorem. 




2 3 
r = R/a 

FIG. 9. The curves shown are the 'minimum' functions, 
9mixx(r)/g, as a function of temperature (such as shown in the 
previous figure) for each value given. These curves are for 
a sc /a = 0.0072. 

Figure [l(^ shows small but significant deviations in the 
Kohn mode from unity for N = 2 • 10 4 and 2 • 10 5 , both 
within DQS-Popov and DQS-G2. That the G2 frequency 
should be lower than the Popov value cannot be simply 
understood in terms of an overall decrease in the inter- 
atomic repulsion, since this would predict a mode closer 
to unity. Rather, the spatial variation of the effective 
interaction leads to a flattening of the effective poten- 
tial, comprised of the trap plus the Hartree potential; 
the looser effective confinement softens all the modes. 
We are not aware of other computational results in which 
the Popov value starts from below unity and rises above, 
before falling near T c . This behavior may be a conse- 
quence of a more rigorous treatment of the chemical po- 
tential, Eq. (|l8|). Alternatively, since the differences in- 
crease with N (specifically, the non-condensate density), 
they may not have been observable with the smaller N 
values studied previously. 

The temperature-dependence of the low-lying excita- 
tion frequencies obtained with the DQSP and DQSG ap- 
proaches is shown in Fig. |Il] for N — 2 • 10 4 and 2 • 10 5 . 
The softening of all the excitation frequencies in the G2 
approximation was found previously by the proponents 
of this theory jll[] (for a 'pancake' geometry) as well as by 
others using a similar perturbative approach to the inter- 
acting Bose gas [Q. However, for a spherically symmet- 
ric trap, the results of Ref. |[l] for 2,000 Rb atoms showed 
only a negligible difference between Popov and G2 ex- 
citation frequencies. The present results show that for 
a spherically symmetric trap and larger atom numbers, 
there can be differences between the Popov and G2 values 
that would be experimentally detectable. These results 
also lead to the question whether for larger atom num- 
bers, a renormalized atom-atom interaction would effect 



frequencies calculated by the methods of Refs. p3| , pT[ , 
which did not assume a static condensate. It should 
also be mentioned that experimentally observed excita- 
tion frequencies with larger numbers of sodium atoms in 
a 'cigar' geometry [Q also exhibited a softening of both 
the quadrupole and dipole excitation frequencies as the 
temperature approaches T c . 
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FIG. 10. Excitation frequencies of the lowest 1 = 1 mode 
in comparison with the Kohn theorem value of unity. Results 
from the Popov (dashed lines) and G2 (solid lines) approxi- 
mations are shown for (a) N = 2 • 10 4 and (b) N = 2 • 10 s . 
All results are for a sc /oo = 0.0072. 
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FIG. 11. Excitation frequencies for the lowest £ — 0, 1, 
and 2 modes for (a) N = 2 • 10 4 and (b) N = 2 • 10 s within 
the Popov (dashed lines) and G2 (solid lines) models, where 
a sc /a = 0.0072. 



V. DISCUSSION AND CONCLUSIONS 

In this work, we have extended finite-temperature 
mean-field calculations for Bose-Einstein condensates 
confined in harmonic traps PJPl[]. A careful derivation 
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of the mean-field equations provides improved definitions 
of the thermodynamic chemical potential and quasipar- 
ticle occupation factors, yielding observables that are 
continuous functions of the particle interactions. The 
numerical techniques employed in the calculations have 
allowed for the investigation of systems with the large 
numbers of atoms relevant to on-going experiments. In 
the process, we have been able to make several crucial 
comparisons between the results of evaluating discrete 
summations over quasiparticle states (which are numeri- 
cally time-consuming) and various local density approxi- 
mations. Furthermore, we have explored the implications 
of a recently proposed gapless theory which takes into ac- 
count pairing correlations. 

The results presented above indicate a significant in- 
adequacy of conventional static mean-field theory for 
computations of excitation frequencies of trapped Bose 
condensates at finite temperatures. For large number 
of atoms and interaction strength, we find appreciable 
deviations of the dipole frequency obtained with either 
the Popov or G2 approximations from expectations of 
the generalized Kohn's theorem. In our computations, 
the condensate is static in the presence of thermal exci- 
tations. The excited dipole mode corresponds approxi- 
mately to out-of-phase motion of the thermal cloud rel- 
ative to the condensate, as observed experimentally 
when the dipole mode of the thermal cloud is excited sep- 
arately. Detailed modelling of such excitation modes has 
been performed only by restrictive parametrization of the 
condensate and thermal cloud in the collisionless Jl4| or 
hydrodynamic |p| regimes. Both of these approaches ad- 
dress the two-fluid nature of these systems, and produce 
dipole modes that satisfy the Kohn theorem exactly. We 
will argue that equilibrium thermal excitations are com- 
puted accurately by the mean-field DQS methods pre- 
sented here. However, any experimental probe of these 
excitations involves perturbative processes that require 
other theoretical methods. 

In principle, mean-field theories that include fluctua- 
tions in the population of excited states |l(],|5l|] ought 
to be equivalent to the two-fluid dynamics in the colli- 
sionless regime. A full second-order perturbation the- 
ory of the interacting Bose gas should yield the cou- 
pled modes of the condensate and thermal clouds as well 
as damping rates. Indeed, employing the approximate 
many-body T-matrix in the calculations (the G2 approx- 
imation described above) yields excitations that have a 
temperature-dependence qualitatively similar to that of 
out-of-phase modes. We hope to explore these issues in 
future work. 
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TABLE I. Comparison of condensate numbers, no, ob- 
tained from Quantum Monte Carlo calculations and from this 
work, with and without atom-atom interactions, and results 
obtained here from discrete Bogoliubov quasiparticle sums 
and discrete Hartree-Fock sums. The error limits for QMC 
are of course purely positive for no = 0. 



T/fuv 

16.667 

16.949 

17.242 

17.544 

17.857 

18.182 

18.519 

18.868 

19.231 

19.608 

19.802 

20.0 

20.202 



QMC 


DQSG 


DQSP 


HFLDA a 


HFLDA b 


N c 


N c 


N c 


N c 


N c 


2265(10) 


2213 


2159 


2216 


2222 


1971(10) 


1936 


1883 


1945 


1935 


1656(15) 


1654 


1599 


1630 


1638 


1374(10) 


1367 


1309 


1323 


1333 


1057(10) 


1072 


1016 


1008 


1022 


741(10) 


782 


726 


686 




440(10) 


501 


448 






180(10) 


247 


205 






21(11) 


140 


57 






0(20) 


71 


21 






0(20) 




15 






0(10) 




12 






0(14) 
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a Holzmann et al. fl7|,[l5| 
b This work, using Eqs. (| 
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